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to 2A^-th Order Initial Value Problem 
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Abstract 



We consider a set of interwoven harmonic oscillators where the accel- 
eration of a given oscillator is determined by the position of its nearest 
' neighbor. We show that this problem of N non-local oscillators with 

periodic boundary conditions leads to a 2A^-th order initial value prob- 
lem. We discuss the numerical solution of this using a non-polynomial 
1-^^ , spline method. A very precise numerical method that minimizes the 

d ' error can be developed, which we test for a few examples of driving 

forces. 

Keywords: Non-local coupled harmonic oscillators, non-polynomial spline; consistency relations; 
end conditions; A^th-order Initial Value Problem. 

> 

(N 

^ ■ 1 Introduction 

(N 

■ Coupled harmonic oscillators are a standard paradigm in many engineering, physical, chemi- 

, cal or biological systems. The basic ingredient of the harmonic oscillator dynamics is the fact 

I that the acceleration of each oscillator is proportional to its position, with a proportionality 

constant given by the negative of the frequency squared. In addition, for non-equilibrium 
situations there are time-dependent driving forces on each oscillator. 
^ I In this paper we are interested in a fundamental modification of this approach. We 

assume that the acceleration of each oscillator is proportional to the position of its neigh- 
boring oscillator. By this, of course, a strongly coupled structure is introduced, and a kind 
of non-local dynamics, since the acceleration of a given oscillator is determined by a posi- 
tion variable elsewhere. Such a dynamics is motivated for very strongly coupled individual 
systems. For example in [H [21 [3] strongly coupled oscillator systems are studied, which 
degenerate to our type of dynamics in the infinite coupling limit. In this limit individual 
local oscillator behavior is infiuenced in a hierarchical way by nearest neighbors. We as- 
sume a driven non-equilibrium situation where each oscillator is also driven by individual 
time-dependent driving forces. 
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We will show that the above dynamics, for N oscillators with periodic boundary con- 
ditions, leads to a high-order initial value problem, indeed of order 2N. Thus rather large 
derivatives will become relevant if there are just a few oscillators coupled in this way. These 
types of initial value problems require effective numerical methods, which depend on the 
number A^. Useful in this context are non-polynomial spline methods. We will describe 
the optimum way to solve this system numerically for a given set of driving forces. In fact, 
we will show that there is an optimum method with minimum error, which is in particular 
useful if highest precision numerical results are required. As main examples, we will deal 
with the cases N = 2 and = 3. 



2 Nonlocal coupling of driven harmonic oscillators 

To illustrate the idea we start with two uncoupled driven harmonic oscillators: 



y\ + u:lyi = gi{t) (2.1) 
2/2 + wiy2 = 52 (i) (2.2) 

LOi is the frequency of oscillator i, and gi{t) is a local driving force on oscillator i. As 
mentioned in the introduction, the basic idea is to consider a modified dynamics where the 
position yi of the first oscillator determines the acceleration force of the second one, and 
vice versa. That is to say, instead of the above trivial (uncoupled) dynamics we consider 
the following coupled dynamics: 

y\+uly2 = gi{t) (2.3) 
m+ojlyi = 52 (t) (2.4) 



Differentiating eq. (12. 4p twice we get 

yf+u:lyf^=g?\t) (2.5) 

where y^^'^ denotes the k-ih. derivative with respect to time t. Eliminating in this equation 
y\ using eq. (|2.3p we get 



+ ul{gi{t) - uly2) = g?\t) (2.6) 
This is equivalent to a 4-th order equation of the form 

+/(t)y2(t) = 5(t) (2.7) 

where 

fit) = -ulu^l (2.8) 

and 

g{t)=gf\t)-ojlgi{t). (2.9) 

Together with the initial position and velocity of the two oscillators, this leads to an initial 
value problem of 4th-order, for which we will describe the optimum numerical method in 
section 3. 
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But let us here first extend the problem, by considering = 3 non-local oscillators, and 
later an arbitrary number A^. The non-locally coupled dynamics for 3 oscillators reads 

y\+ujly2 = gi{t) (2.10) 
m+ujiys = g2{t) (2.11) 
m+ujlvi = g3{t) (2.12) 

Differentiating eq. (I2.12|) twice we get 

yr+a.i,f)=,f(i). (2.13) 
Eliminating in this equation y\ using eq. (j2.10p we get 

yi'^ + co|{g^{t) - coly2) = g^^\t) (2.14) 



or 



yi^^ - ojluly2 = gf\t) - ulgi{t). (2.15) 



Differentiating twice this leads to 

yf-u:lulyf^=gf\t)-ulgf\t) (2.16) 
and using (|2.1ip to eliminate we get 

yf + ujlululy^i = gf{t) - ujlgf\t) + wiw?52(t). (2.17) 
Apparently, this can be written in the form 

yf +f{t)y:,{t)=g{t) (2.18) 

where 

f{t) = (2.19) 

and 

9{t) = gf (t) - c^isf ^ (t) + ^liolg2{t). (2.20) 

It is obvious how to generalize this problem to A^ oscillators. In this case one obtains a 
2A^-th order initial value problem of the form 

yS''^ + fit)yN = git) (2.21) 

where 

N 

/(t) = (-1)^+1 (2.22) 

i=l 

and g{t) is a sum of various derivatives of local driving forces gi{t) weighted with frequencies. 

The optimum numerical way to solve these high-order initial value problems depends 
on A'^ in a nontrivial way. In the following, we allow for general time-dependent functions 
f{t) and deal in detail with the cases N = 2 and A^ = 3. 
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3 N=2: Optimized numerical solution of the 4th order initial 
value problem 

3.1 The initial value problem 

The case of two oscillators {N = 2) leads to the following fourth order initial value problem 

via) =uo, y^'Ka) = ui, } (3.23) 

y(2)(a) = U2, 

where the Ui{i = 0,1,2,3) are finite real constants while the functions f{t) and g{t) are 
continuous on [a,b]. To simplify the notation, we have written y{t) instead of 2/2 (t). 

Of course, from an engineering point we are interested in four initial values given by 
initial position and velocity of oscillator 1 and 2. These are related to the constants Uj by 

y2{a) = uo (3.24) 

2/2(0) = ui (3.25) 

32 (a) - U2 , . 

yi(«) = 2 — ("^-^S) 

. , ^ 32 (g) -^^3 

2/i(«) = 2 (3-27) 

^2 

We have reduced the 4-dimensional problem of the space-space structure of the two coupled 
nonlocal oscillators to a 1-dimensional initial value problem of 4-th order, for which we can 
apply very precise numerical methods, as described in the following subsection. While for 
our nonlocal oscillator example f{t) = const = — w^wl, the numerical method developed 
in the following is applicable general time-dependent functions /(t), as long as they are 
continuous. 



3.2 Nonpolynomial Spline Method 

To develop the spline approximation to the problem (j3.23p . the interval [a, h] is divided into 
n equal subintervals, using the grid points U = a + ih ; i = 0, 1, . . . , n, where h = {b — a)/n. 
Consider the following restriction Si of the solution to each subinterval [tj, tj+i], i = 
0,l,...,n-l, 

Si{t) = aiCOSLo{t - ti) + bismu}{t - U) + Ci{t - tif + di{t - tj)^ -|- ei{t - U) + pi 

(3.28) 

Let 

Vi = Siiti) Mi = Si^\ti), 

Ni = 5f)(t.), 

Following [3] and postulating that at the end points of the intervals the 1st and 3rd deriva- 
tives are continuous, one derives the following consistency relation between the values of 



0,1, 



(3.29) 
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splines and their fourth order derivatives at border points: 

= [Vi-A - %-3 + 6yi_2 - 4yi_i + yt] ; 

z = 4,5,...,n, 



(3.30) 



where 



a 



1 



1 



e^sin^ ^Ssin^ 



1 

+ 777 



2(1 + cos ( 
6*3 sin 6* 



(cos 9-2) 4 



3^ 



and 



7 



-2(1 + 2 cos 61) (1-4 cos 6*) 6 
^3 sine ^ 30 sine ^ ' 



Here 6 = uh is an arbitrary parameter. The relation (|3.30p forms a system of n — 3 linear 
equations in the n unknowns 

{Ui, i = 1,2, ...,n), while Ni is taken from IVP (j3.23p to be equal to —fiin + gi, 
i = 0,1, . . . ,n. 

Following [5], three equations (end conditions) are determined to find the complete solution 
of yiS appearing in eq. (I3.30[) . as given below: 



N0 + N4 



1 



-220 , 40 40, (1) 4,4 (4) 

-^yo + 40yi - 20y2 + y 2/3 - y %r " ^h^o 



18336 
575 



-yi 



22992 
575 



-y2 



4656 
'575 



■y-i + 



2736 
115 



(1) 



575 



575 



(3.31) 



(3.32) 



and 



N2 + NQ 



1 



8157 11424 3267 978, (1) 8958 (2) 
y2 - ^n^ys + -^yi + -r^hy'^ + -— h y'o 



865 



865 



865 



173 



865 



+ 5684 3 (3) 
+ 865 " ^0 



(3.33) 



The local truncation errors associated with the linear equations (j3.3ip — ()3.33p and (j3.30p 
are calculated as 



-fh^y(^Hh)+0{h'^), 

-^/.V^)(t2) + 0(/.^), 
.l^h^yi^)it,)+0{h'). 



i = 1 
i = 2 
i = 3 



i(-l + 24a + 6/3)/i6y(6)(t.) + 0{h'^), i = 4, 5, 



and 



|r|| = ch'^Ri = 0{h^ 



Ri = max ly'-^Xt) 

te[a, b] 



(3.34) 



(3.35) 



where c is a constant which depends only on the values of a and (3 and is independent of 
h. 
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Let us mention that the solution obtained using the system of hnear equations (j3.3ip — 
(j3.33p and (|3.30p is second order convergent. But if a, /3 and 7 are taken such that 
^ ~ ~750' ~ MJ' ~ ^ then the order of the truncation error in eq. ()3.30p is 
0{h^^) and the order of convergence can then be improved up to sixth order, using this 
method of improved order of end conditions. 



The improved end conditions with truncation error of order 0{h ) are 



1 

7? 



7? 



-68397280 
18729 i'O 



^ , 843268 ^ , 330342 ^ _ 16892 ^ , ^ 
^^0^ 2081 2081 2081 ^JV4 

13366080 ., _ 7408800,, , 14781760,, _ 10427200 , 743680 1,2„(2) 

2081 2081 y2 -I- 18729 i/3 5243 "-^/o 2081 ^0 

AT _ 156090207332 at _ 40456201386 at _ 600708692 at 1 /V 
158360705 ^^2 158360705 ^^3 158360705 "T 

180155114496 _ 340726283352 210168798336 „ _ 49597629480 „, , 69181575120 l„(1) 

31672141 31672141 ^2 "I" 31672141 ^3 31672141 ^4 "I" 31672141 "-^O 

I 42396452784 ;,2„(2) , 7557647328 A,3„(3)l /o o'j^ 

^ 31672141 " yO "T 31672141 " ^0 J 



(3.36) 



and 



1 

7? 



N2 



85514900495708 
1252977040745 



A^3 + 



3759590586966 
1252977040745 



A^4 



7418340285788 
1252977040745 



43463161469952 
250595408149 ^2 



94491207986112 , 68699611790208 „, 
250595408149 ^3 250595408149 ^4 



17671565274048 „ 
250595408149 ^5 



10106680227840 



9581784601536 f,2 (2) , 2621304758016 l 3 „, (3) 
250595408149 "i^O 250595408149 i^O 250595408149 ^ ^0 



(3.38) 



The truncation errors of the corresponding equations are 

-0.3034/ii0y(i0)(ti) + O(/iii), 
-1.4034^i%(i0)(t2) + O(/i^^), 
-1.0163/ii%(i0)(t3) + O(/iii), 
^{-17 + 5376a + 8i/3)h^^y(^°\ti) + 0{h 



11^ 



1, 
2, 
3, 

4,5, 



(3.39) 



and 



T|| = ch''^R2 = 0{h 



10\ 



R2 



max \y 

tG[a, b] 



(10) 



(3.40) 



where c is a constant which depends only on the values of a and /? and is independent of 
h. 

To illustrate the powerfulness of the method, two analytically solvable examples are dis- 
cussed in the following: 



3.3 Examples 



Example 1 
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Consider the following IVP 



y 



y 

(3) 



yit) = 4cos(t), t e [-1, 1], 

-1) = -2sm(l), 

-1) = 2cos{l) + sin{l), 

-1) = -2cos(l) + 2sm(l) 

-1) = -2cos(l) - 3sm(l). 



(3.41) 



The analytic solution of the above problem is 

y{t) = (l-t) sin{t) . 

The observed maximum errors (in absolute values) associated with yi, for the problem 
(|3.4ip . corresponding to the different values of a, (3 and 7, are summarized in Table 1. It is 
confirmed from Table 1 that if h is reduced by factor 1/2, then \\E\\ is reduced by a factor 
1/4, which indicates that the method gives second-order results. 



Table 1: Maximum absolute errors for problem (j3.4ip in y^. 



n 


a = 0, /3 = 


a = 1/2, /? = 1/2 


a = 1/6, /? = 1/6 




7 = 1 


7 = -l 


7 = 1/3 


6 


6.74 X 10"^ 


3.6 X 10° 


1.73 X 10° 


12 


5.77 X 10-2 


7.3 X 10-1 


2.22 X 10-1 


24 


3.3 X 10-3 


4.5 X 10-2 


1.3 X 10-2 


48 


1.48 X 10"^ 


2.1 X 10-3 


5.93 X 10-"^ 
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The observed maximum errors (absolute values) associated with yi for the problem 
(j3.4ip . corresponding to the use of improved end conditions, are summarized in Table 2. 
A significant improvement of precision is obtained. 



Table 2: Maximum absolute errors for problem (j3.4ip in yj. 



n 


\y{U) - yi\ 


6 


1.7 X 10-3 


12 


1.17 X 10-^ 


24 


7.19 X 10"^ 


48 


7.72 X 10"^^ 



Example 2 



Consider the following IVP 



+ty{t) 
y(0) = 0, 
y(2)(0) =0, 



+ 7t + t3), 0<t<l 



y«(0) 
y(3)(0) = 



-- 1, 
-3. 



(3.42) 



The corresponding analytic solution is now 

y{t) = t{l -t)e' . 

The observed maximum errors for different values of a, /3 and 7 are summarized in Table 
3. It is confirmed from Table 3 that if h is reduced by factor 1/2, then \\E\\ is reduced by 
a factor 1/4, which indicates that the method gives second-order results. 
The observed maximum errors using improved improved end conditions are summarized in 
Table 4. 



4 N=3: Optimized numerical solution of the 6th order initial 
value problem 

4.1 Non-polynomial spline solution for N=3 

A similar optimized numerical method can be developed for the case = 3. In this case 
the IVP reads 
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Table 3: Maximum absolute errors for problem (j3.42p in yj. 



n 


a = 0, /3 = 


a = 1/2, (3 = 1/2 


a = 1/6, (3 = 1/6 




7 = 1 


7 = -l 


7 = 1/3 


6 


1.14 X 10"^ 


2.31 X 10^2 


6.86 X 10-2 


12 


1.14 X 10~2 


1.55 X 10^2 


2.4 X 10-3 


24 


1.4 X 10-3 


4.8 X 10^3 


6.40 X IQ-'' 


48 


2.18 X 10"^ 


1.3 X 10^3 


2.87 X 10-"^ 



Table 4: Maximum absolute errors for problem (j3.42p in yj. 



n 


\y{k) -yi\ 


6 


2.53 X 10-^ 


12 


1.53 X 10""^ 


24 


1.06 X 10-9 


48 


1.09 X 10-^° 
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(4.43) 



y^'\t) + f{t)y{t) = g{t), t G [a, b], 
y{a) = no, y(i)(a) = ni, 

y{2)(a)=ti2, y(3)(a)=^i3, 

y(4)(a) = n4, y(^)(a) = ^5, 

Again the interval [a, 6] is divided into n equal subintervals, using the grid points ti = a + ih 
(i = 0, 1, n), where /i = {h — a)/n. 

Again we consider the restrictions Si of the solution to each subinterval [tj, tj+i], i = 
0,1,... ,n- 1, 

Si{t) = Qi cos uj {t - ti) + bi sin uj{t - ti) + Ci{t - ti)^ + di{t - ti)'^ 

+ei{t - tif + qi{t - tif + ri{t - U) + Vi. (4.44) 



0,1,..., n. (4.45) 



and define 

= Si{ti) M, = Sf\u), 

N, = Sf^\u), U = sf\ti), 

We denote by y{t) the exact solution of the IVP (|4.43p and yi is the approximation to y{ti), 
obtained by the spline S{ti). Prom continuity of the first, third and fifth derivatives at the 
border points, i.e. sl^i{ti) = sl^\ti), fi= 1,3 and 5, one gets 

L6r f e-smd 1 I 1 

*-6 V e** sine eP' 
4_/,6r f ^ — 2(cos 9+2) cosg-l _ cos 0-13 

,,6r / (8cose+7) 15 I 4cos6>+5 (52 cos 0-67) \ 

sinO 12Of?sin0 / 

_L/,6r /20 2(6cos6>+4) 2(3cos6>+l) 33 cos 0-13 ' 
-rri. 1^1-3 \gU e^sine Se-i sine 3O6»sm0 

,6r / (8cose+7) 15 I 4cos6>+5 (52 cos 6>-67) \ ' (4.46) 

sin 8 120s sin 6* / 

4_^6r / 6 _ 2(cos6l+2) cos 0-1 _ cos 0-13 A 
-I-"- -^j-1 l^gB- e5sin6» 363 sine GOesine^ 

-Lh^T ■ ( e-sine _ 1 , 1 A 

« V e^sine 69^ sin 9 12^ sin ^ / 

= [yi-6 - 6yj-5 + 152/j-4 - 20?/i_3 + 15yi-2 - Qyi-i + Vi] ; 
i = 6,7, . . . ,n, 

which can further be written as 

{ah^Li.e + /J/i'^ii-s + lh^Li^4 + 5/i^ii-3 + lh<^L,^2 + fih^ U^i + a/i^L^) 
= [yi-Q - 6yi-5 + 15yi_4 - 20yi_3 + 15yi-2 - ^Vi-i + Vi] ] 

i = 6,7,...,n, (4.47) 



where 



a 



sine 1 1 



6*6 sin 61 66'3sin6l 1261 sin ( 



_ 6 2(cos6l + 2) cos6'-l cos6'-13 

^ ~ ^¥ e^sine sin 9 609 sin 9 

(8cos6' + 7) 15 4cos6' + 5 (52cos6'-67) 

^ ~ ' 615 sin 61 F"^ 66*3 sin 6* 1206* sin6l 
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and 



20 



2(6 cos (9 + 4) 2(3cos6' + l) 33cos6'-13 



9^ sine 3(93 sin 6* 30(9 sin 6^ 

Here 9 = uh. The relation (I4.47P forms a system of n — 5 linear equations in the n unknowns 
{Ui, i = 1,2, n), while Lj is taken from IVP (|3.23p to be equal to —fiUi + Qi, 
i = 0,1, . . . ,n. 

Following [6], five equations (end conditions) are determined to find the complete solution 
of UiS appearing in eq. (I4.47[) . as given below: 



Ln + La 



1 



1 

hp 



2905 



12 



-yo - 336yi + 126^2 



112 21 , (i) 



+42/.2yf 



797790 
21983 



1660890 



(4.48) 



■yi 



87150 



1 

hp 



21983 
605725 
22267 
675200 

+ 



21983 
283500 
21983 



y2 + 



1299060 



21983 



-ys 



523110 
21983 



hyil' + 



108239440 
1803627 



-ys + 



172620 
21983 
1103910 
22267 



2„(2) 



h'y'o 



■yA 

40167 



21983 
446800 5949805 



-yA 



22267 



■ys + 



1803627 



(4.49) 



ye 



(1) 700180 ^ (2) 851440 ^3 (3) 



L3 + L7 — — 7- 



1 



85887 " 66801 

670672000 44149995 

ys + T77::TTT7;^y4 



42346017 
12961750, (1) 

hyo 



2016477 
11492010 ^4 (4) 

1568371 



1568371 

25078370 ,0 (2) 

^ yo 



200403 

23862240 122902615 
ys + -n^TTTTTTTTrye 



(4.50) 



1568371 



1568371 " 42346017 
77684300 
4705113 ^° 



(4.51) 



and 



L4 + L« 



1 

7^ 



49567095 
12837314^"^ ' 
59244435 

+■ 



34289280 



6418657 



■ys + 



19011465 



12837314 



ye + 



2182545 
916951 



■hy^ 



(1) 



6418657 



h'y'^' + 



107795790 
6418657 



/.^yf + 



115282605 



4„(4) 



6418657 



-h'^y'o 



65492262 . (5) 



6418657 

The local truncation errors associated with the linear equations (j4.48p 
are calculated, as 

' -4.75/i8y(8)(ti) + 0(/i9), i = 1, 

-5.0467/i8y(8)(t2) + 0{h'^), i = 2, 

-5.9909/i8y(8)(t3) + 0(/i9), i = 3, 

-12.3201/i8y(s)(t4) + 0(/i9), i = 4, 

-23.7869/i8y(8)(t5) + 0{h^), i = 5, 
(-1 + 2a + 2/3 + 27 + 6)h^y'^^\ti) 

+i(-l + 36a + 16/3 + 47)/i8y(8)(ti) + 0{h^), i = 6,7,.. 



(4.52) 

(|i32D and dOT]) 



(4.53) 



, n. 
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To make the truncation errors of the system (|4.47p of order h^, a, /3, 7 and 6 are taken such 



that a + /3 + 7 + | = i and then 



jr|| = ch^Rs = 0{h^ 



where c is a constant. 



i?3= max 

ie[a, h] 



(4.54) 



The solution obtained using the system of Hnear equations (j4.48p — (j4.52p and (j4.47p in 
general is second order convergent. Again, however, the order of accuracy of the method 
can be improved significantly to . The local truncation error of the system (j4.47p can be 
expressed in the following form 



1 + 36a + 16/3 + ^-i)h^y^^\ti 



(4.55) 



(-1 + 2a + 2^ + 27 + 5)h^y^'^\ti) + ^y-i 
+ 230 (-7 + 1620a + 320/3 + 20j)h^^y'^^^\tij 
+ 7^(-16 + 15309a + 1344/3 + 21 j)h^^y^^^\ti) 
+mo (-13 + 39366a + 1536/3 + 6^)h^^y^^^'^ (U) 
+ 159667200 (-6^1 + 5196312a + 90112/3 + 88-f)h^^y^^^\t 

i = 6,7, ... ,n, 

Thus, the order of the truncation error ii can be improved to be of order h^^ and corre- 
spondingly the order of method can be improved up to /i^ , if a = 3054Q , /3 = 41/5040, 7 = 
2189/10080, 6 = For other choices of the parameters (not listed here), one can make 
the method to be of order h?,h^,h^, respectively. Results corresponding to the order h?, 
h^, and are described in the following section. 



4.2 Test of the method with analytically solvable examples 



Example 3 

Consider the IVP 



y(6)(t) 


-y{t) = 


-6e*, < 


t < 1 


y(o) 


= 1, 


y{l)(0) 


= 0, 


y(2)(o) 


= -1, 


y(3)(o) 


= -2, 


yW(o) 


= -3, 


y(5)(0) 


= -4. 



(4.56) 



The analytic solution of this IVP ()4.56p is 

y{t) = (1 - t) e* . 

The observed maximum errors are summarized in Table 5. It is confirmed from Table 5 
that if h is reduced by factor 1/2, then \\E\\ is reduced by a factor 1/4, which indicates that 
the method gives second-order results. 

The observed maximum errors (in absolute values) associated with yi, for the problem 
()4.56p . corresponding to different orders of method are summarized in Table 6. 
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Table 5: Maximum absolute errors for problem (j4.56p in yj. 



n 


a = 1/120, P = 15/120 
7 = 1/4, <5 = 28/120 


" 720' ^ 36 

7 = 219 (J = 240 
' 720' 720 


" 5040' ^ 504 

1250 X _ 2418 
' 5040' " 5040 


8 


7.98 X 10-"^ 


9.13 X 10"^ 


9.51 X 10-^ 


16 


7.50 X 10-^ 


9.64 X 10"^ 


1.03 X 10-^ 


32 


5.45 X 10-6 


1.02 X 10-5 


1.18 X 10-5 


64 


1.28 X 10"^ 


9.42 X lO-'^ 


1.37 X 10-6 



Table 6: Maximum absolute errors for problem (j4.56p in 



n 


0{h') 


0(/i6) 


0{h^) 


8 


4.04 X 10-5 


2.07 X 10-^ 


2.13 X IQ-i 


16 


1.10 X 10-6 


8.99 X 10-9 


4.80 X lO-'^ 
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Example 4 

The second example is 

+y{t) = 6(2tcos(t) + 5sin(i)), -1 < t < 1 
yi-1) =0, y^^\-l) =2sm(l), 

y(2)(_i) = -4cos(l) - 2sm(l), = 6cos(l) - 6sm(l), 

y(4)(_i) = 8cos(l) + 12sm(l), 2/^^^(-l) = -20cos(l) + lOsm(l). 

with analytic solution 

y{t) = {t^ - 1) sin(t) . 

The observed maximum errors (in absolute values) associated with yi, for the system (|4.57|) . 
corresponding to the different values of a, /3 7 and 5 are summarized in Table 7. Again 
it is confirmed from Table 7 that if h is reduced by factor 1/2, then \\E\\ is reduced by a 
factor 1/4, which indicates that the method gives second-order results. 
The observed maximum errors corresponding to the different orders of the method are 



Table 7: Maximum absolute errors for problem (I4.57P in yi. 



n 


a = 1/120, (3 = 15/120 
7 = 1/4, (5 = 28/120 


~ 720 ' ~ 36 

219 5-240 
' 720' " 720 


'-^ ~ 5040 ' ^ ~ 504 

1250 X _ 2418 
1 — 5040' " ~ 5040 


16 


7.35 X 10-2 


9.64 X 10-2 


1.03 X 10-1 


32 


1.01 X 10-2 


1.62 X 10-2 


1.82 X 10-2 


64 


4.51 X 10""^ 


2.0 X 10-3 


2.5 X 10-3 


128 


1.98 X 10-^ 


1.79 X lO-'' 


3.05 X 10-"^ 



summarized in Table 8. 



5 Conclusion and Outlook 

In this paper we started from a set of N nonlocal coupled harmonic oscillators, each driven 
by a driving force. We showed that this leads to an 2A^-th order initial value problem 
(IVP) in a single variable. In a sense this gives 'physical meaning' to high-order IVP in 
one variable, which so far have mainly been looked at without any physical interpretation. 
Engineering applications include strongly coupeld oscillator problems where the state of a 
local oscillator is strongly influenced by the position of the nearest neighbor. 
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Table 8: Maximum absolute errors for problem ()4.57p in m. 



n 


0{h') 


0{h^) 


0{h^) 


8 


2.31 X 10^2 


2.87 X 10"^ 


2.98 X IQ-^ 


16 


8.6 X 10-3 


7.98 X 10-5 


9.93 X 10-8 



By implementing improved end conditions, a very precise numerical method could be 
developed to solve this system numerically. In fact, we believe it is one of the most pre- 
cise methods known in the field. Apparently our results are relevant to find very precise 
numerical solution schemes for higher-dimensional differential equations. We showed that 
a transformation of an A^-dimensional 2nd order differential equation to a 1-dimensional 
differential equation of order 2N can be highly advantageous from a numerical point of 
view. One can implement improved end conditions that allow for a significant reduction of 
the error. After the 2A^-th order IVP has been solved very precisely, the solution can be 
translated back into the original physical setting of nonlocal oscillators. 

While we have explicitly worked out the cases N = 2 and A^ = 3, in principle our 
method can be extended to higher values of A^, though the complexity of the formulas used 
to minimize the truncation error increases rapidly. 
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